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We present a new approach to simulating Hamiltonian dynamics based on implementing linear 
combinations of unitary operations rather than products of unitary operations. The resulting algo- 
rithm has superior performance to existing simulation algorithms based on product formulas and, 
most notably, scales better with the simulation error than any known Hamiltonian simulation tech- 
nique. Our main tool is a general method to nearly deterministically implement linear combinations 
of nearby unitary operations, which we show is optimal among a large class of methods. 



Simulating the time evolution of quantum systems is a major potential application of quantum computers. While 
quantum simulation is apparently intractable using classical computers, quantum computers are naturally suited to 
this task. Even before a fault-tolerant quantum computer is built, quantum simulation techniques can be used to prove 
equivalence between Hamiltonian-based models of quantum computing (such as adiabatic quantum computing [i] and 
continuous-time quantum walks [2]) and to develop novel quantum algorithms [3-7]. 

In recent years there has been considerable interest in optimizing quantum simulation algorithms. The original 
approach to quantum simulation, based on product formulas, was proposed by Lloyd for time-independent local 
Hamiltonians [8] . This work was later generalized to give efficient simulations of sparse time-independent Hamiltonians 
that need not have tensor-product structure [4, 9]. Further refinements of these schemes have yielded improved 
performance [lO-l.'-!] and the techniques have been extended to cover time-dependent Hamiltonians [14, 15]. Recently 
a new paradigm has been proposed that uses quantum walks rather than product formulas [IG, 17]. This approach is 
superior to the product formula approach for simulating sparse time- independent Hamiltonians with constant accuracy 
(and in addition, can be applied to non-sparse Hamiltonians), whereas the product formula approach is superior for 
generating highly accurate simulations of sparse Hamiltonians. 

The performance of simulation algorithms based on product formulas is limited by the fact that high-order ap- 
proximations are needed to optimize the algorithmic complexity. The best known high-order product formulas, the 
Lie-Trotter-Suzuki formulas, approximate the time evolution using a product of unitary operations whose length 
scales exponentially with the order of the formula [IS]. In contrast, classical methods known as multi-product for- 
mulas require a sum of only polynomially many unitary operations to achieve the same accuracy [ 9] (although of 
course the overall cost of classical simulations based on multi-product formulas remains exponential in the number of 
qubits used to represent the Hilbert space). However, these methods cannot be directly implemented on a quantum 
computer because unitary operations are not closed under addition. 

Our work addresses this by presenting a non-deterministic algorithm that can be used to perform linear combinations 
of unitary operators on quantum computers. We achieve high success probabilities provided the operators being 
combined are near each other. We apply this tool to quantum simulation and thereby improve upon existing quantum 
algorithms for simulating Hamiltonian dynamics. Our main result is as follows. 

Theorem 1. Let the system Hamiltonian be H = X^JLi where each Hj G is Hermitian and satisfies 

^ ^ for a given constant h. Then the Hamiltonian evolution e^*^* can be simulated on a quantum computer 
with failure probability and error at most e as a product of linear combinations of unitary operators. In the limit of 
large m, ht, 1/e, this simulation uses 



elementary operations and exponentials of the Hj s. 

Although we have not specified the method used to simulate the exponential of each Hj, there are well-known 
techniques to simulate simple Hamiltonians. In particular, if Hj is 1-sparse (i.e., has at most one non-zero matrix 
element in each row and column), then it can be simulated using 0(1) elementary operations [4, !)], so (1) gives an 
upper bound on the complexity of simulating sparse Hamiltonians. 

Our simulation is superior to the previous best known simulation algorithms based on product formulas. Previous 
methods have scaling of the same form, but with the coefficient 1.6 replaced by 2.54 [11, Theorem 1] or 2.06 [20, 
Theorem 1]. Also note that Theorem 1 of [12] gives a similar scaling as in [20], except the term in the exponential 
depends on the second-largest \\Hj\\ rather than h. 



I. INTRODUCTION 
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Perhaps more significant than the quantitative improvement to the complexity of Hamiltonian simulation is that 
our approach demonstrates a new class of simulation protocols going beyond the Lie-Trotter~Suzuki paradigm, the 
approach used in most previous simulation algorithms. It remains unknown how efficiently one can perform quantum 
simulation as a function of the allowed error e, and we hope our work will lead to a better understanding of this 
question. 

The remainder of this article is organized as follows. In Section II, we provide a general method for implementing 
linear combinations of unitary operators using quantum computers and lower bound its success probability. This 
method is optimal among a large class of such protocols, as shown in the appendix. In Section III, we provide a brief 
review of Lie-Trotter-Suzuki formulas and multi-product formulas and then show how to implement multi-product 
formulas on quantum computers. Error bounds and overall success probabilities of our simulations are derived in 
Section IV. We then bound the number of quantum operations used in our simulation in Lemma 12, from which 
Theorem 1 follows. We conclude in Section V with a summary of our results and a discussion of directions for future 
work. 



II. ADDING AND SUBTRACTING UNITARY OPERATIONS USING QUANTUM COMPUTERS 

In this section we describe basic protocols for implementing linear combinations of unitary operations. Lemma 2 
shows that a quantum computer can nearly deterministically perform a weighted average of two nearby unitary 
operators. (Our approach is reminiscent of a technique for implementing fractional quantum queries using discrete 
queries [ ].) We build upon Lemma 2 in Theorem 3, showing that a quantum computer can non-deterministically 
implement an arbitrary linear combination of a set of unitary operators. 

Lemma 2. Let Ua, C4 € he unitary operations and let A = \\Ua — Ub\\- Then for any k > 0, there exists a 

quantum algorithm that can implement an operator proportional to KUa + Ui, with failure probability at most A^k/(k + 
1)2 < in/iii+iy. 

Proof. Let 

:= r 7^ I • (2) 




Our protocol for implementing the weighted average of Ua and Ub works as follows (see Figure 1). First, we perform 
V„ on an ancilla qubit. Second, we perform a zero-controlled Ua gate and a controlled Ui, gate on the state \tfj) using 
the ancilla as the control. Finally, we apply to the ancilla qubit and measure it in the computational basis. This 
protocol performs the following transformations: 
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If the first qubit is measured and a result of is observed, then this protocol performs \4>) i-> {nUa + Ub)\ip) (up to 
normalization). If the measurement yields 1 then the algorithm fails. The probability of this failure, P+, is 

p ^ \\Ub-Ua\\'^ ^ 
+ - (k+1)2 (At +1)2- ^' 

Since A < 2, this is at most , '^"iso . □ 

By substituting Ub —Ub, Lemma 2 also shows that unitary operations can be subtracted. (Alternatively, replacing 
with K in Figure 1 also simulates Ua — Ub-) Similarly, we could make the weights of each unitary complex by 
multiplying Uq and Ui by phases, although we will not need to make use of this freedom. 

General linear combinations of unitary operators can be performed by iteratively applying Lemma 2. The following 
theorem gives constructs such a simulation and provides bounds on the probabilities of failure. 
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FIG. 1: Quantum circuit for non-deterministically performing an operator proportional to nUa + Ut given a measurement 
outcome of zero. 



Theorem 3. Let V Y.\tl CgUg for k > 1 where Cq ^ 0, \\Ug\\ = 1, and max,^,/ \\Ug - Ug'\\ < A. Let k := 
(Eg: c„>o'~^<])/(J2q: c„<o l^qD- Then there exists a quantum algorithm that implements an operator proportional to 

V with probability of failure + P_ with 

fcA2 

P+ < —, (5) 
4k , , 

where is the probability of failing to add some pair of operators and P_ is the probability of failing to perform the 
subtraction. 

Proof. Let 

We implement V oc — i? using circuits that non-deterministically implement operators proportional to A and B. 
We recursively perform these sums by using the circuits given in Lemma 2, except that we defer measurement of the 
output until the algorithm is complete. We then implement kA ~ B hy using our addition circuit, with Ua taken to 
be the circuit that implements A and Ub taken to be the circuit that implements —B. Then we measure each of the 
control qubits for the addition steps. Finally, if all the measurement results are zero (indicating success), we measure 
the control qubit for the subtraction step. If that qubit is zero, then we know from prior analysis that the non-unitary 
operation V is implemented successfully. 

The operators A and B are implemented by recursively adding terms. For example, the sum Ui + U2 + + U4 is 
implemented as {{{Ui + U2) + + U4). Implementing the sums in A and B requires k — \ addition operations, so 

V can be implemented using k — \ addition operations and one subtraction operation (assuming that all the control 
qubits are measured to be zero at the end of the protocol) . 

According to Lemma 2, the probability of failing to implement V ^ given that we successfully implement A and _B, 

is 

The probability of failing to perform the fc — 1 sums needed to construct A and B obeys 

^-^('-'>(^^"- <'°' 

where the last step follows from the fact that we can take k > 1 without loss of generality. □ 

These results show that we can non-deterministically implement linear combintations of unitary operators with high 
probability provided that k 3> 1 and A <C 1. As we will see shortly, this situation can naturally occur in quantum 
simulation problems. 

Note that it is not possible to increase the success probability of the algorithm by replacing the single-qubit unitary 
V„ with a different unitary, even if that unitary is allowed to act on all of the ancilla qubits simultaneously. We 
present this argument in Appendix A. 
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III. IMPLEMENTING MULTI-PRODUCT FORMULAS ON QUANTUM COMPUTERS 

In this section we present a new approach to quantum simulation: we approximate the time evolution using a 
sequence of non-unitary operators that are each a linear combination of product formulas. Such sums of product 
formulas are known in the numerical analysis community as multi-product formulas, and can be more efficient than 
product formulas for classical computations [19, 22]. We show how to implement multi-product formulas using 
quantum computers by leveraging our method for non-deterministically performing linear combinations of unitary 
operators. 



A. Review of Lie— Trotter— Suzuki and Multi-Product Formulas 



Product formula approximations can be used to accurately approximate an operator exponential as a product of 
operator exponentials that can be easily implemented. Apart from their high degree of accuracy, these approximations 
are useful because they approximate a unitary operator with a sequence of unitary operators, making them ideally 
suited for quantum computing applications. 

The most accurate known product formula approximations are the Lie-Trotter-Suzuki formulas, which approximate 
g-iHt j^Qj. _ Hj as a product of the form 

fc=i 

These formulas are recursively defined for any integer x > by [ ] 

m 1 
j=l j=m 

S^{t) = {S^^,{s^^,t)f 5x_i([l - 4s^_i]i) {S^^,{sx-ityf , (11) 

where Sp = (4 — 4^/'^^'+-'^^)^^ for any integer p > 0. This choice of Sp is made to ensure that the Taylor series of 
matches that of e~^^* to 0{t^^'^^). Consequently, the approximation can be made arbitrarily accurate for suitably 
large values of x ^nd small values of t. 

The advantage of these formulas is clear: they are highly accurate and approximate U{t) :— e~'^* by a sequence of 
unitary operations, which can be directly implemented using a quantum computer. The primary disadvantage is that 
they require 0(5'^) exponentials to construct an 0{t'^''~^^) approximation. This scaling leads to quantum simulation 
algorithms with complexity A product formula requiring significantly fewer exponentials could result 

in a substantial performance improvement over existing product formula-based quantum simulation algorithms. 

In the context of classical simulation, multi-product approximations were introduced to address these problems [19, 
22]. Multi-product formulas generalize the approximation-building procedure used to construct to allow sums 
of product formulas. The resulting formulas are simpler since it is easier to construct a Taylor series by adding 
polynomials than by multiplication alone. Specifically, multi-product formulas only need 0{k^) exponentials to 
construct an approximation of U{t) to 0(t^'^+^). 

We consider multi-product formulas of the form 

fe+i 

M,,^{t)^Y.^,S^(^/i,Y^, (12) 

9=1 

where the formula is accurate to 0{t'^^''~^^^'^^). Here ii,. . . , £k+i are distinct natural numbers and Ci, . . . , Ck+i G K 
satisfy X]g=i — 1- Explicit expressions for the coefficients Cq are known for the case x = 1 [22]: 

n (13) 

j = {l,...,fe+l}\q 1 3 

We will show later that the same expressions for Cq also apply for x > 1. 

For classical simulations, the most numerically efficient formulas correspond to iq — q and x = 1. The simplest 
example of a multi-product formula is the Richardson extrapolation formula for Si [ ] : 

uit)='-^^m^^oit% (14) 
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Even though the above expression is non-unitary, it is very close to a unitary operator. In fact, there exists a unitary 
operator within distance 0{t^^) of the muhi-product formula. In general, Blanes, Casas, and Ros show that if a multi- 
product formula is accurate to 0(t^('^+'^)+^) then it is unitary to 0(t^(*^+^)+^) [19]; therefore such formulas are 
practically indistinguishable from unitary operations in many applications [19, 22]. 

The principal drawback of these formulas is that they are less numerically stable than Lie-Trotter-Suzuki formulas. 
Because they involve sums that nearly perfectly cancel, substantial roundoff errors can occur in their computation. 
Such errors can be mitigated by using high numerical precision or by summing the multi-product formula in a way 
that minimizes roundoff error. 

An additional drawback is that linear combinations of unitary operators are not natural to implement on a quantum 
computer. Furthermore, our previous discussion alludes to a sign problem for the integrators: the more terms of the 
multi-product formula that have negative coefficients, the lower the success probability of the implementation of 
Theorem 3. This sign problem cannot be resolved completely because, as shown by Sheng [24], it is impossible to 
construct a high-order multi-product formula of the form in (12) without using negative Cq. Nevertheless, we show 
that this problem is not fatal and that multi-product formulas can be used to surpass what is possible with previous 
simulations based on product formulas. 



B. Implementing Multi-Product Formulas Using Quantum Computers 

We now discuss how to implement multi- product formulas using quantum computers. The main obstacle is that 
the multi-product formulas most commonly used in classical algorithms have a value of k (as defined in Theorem 3) 
that approaches 1 exponentially quickly as k increases. Thus the probability of successfully implementing such multi- 
product formulas using Theorem 3 is exponentially small. 

Instead, we seek a multi-product formula M^^^, with a large value of k, such that ||Mfc^^(A) — U{X)\\ G 0(A'^('^+^)+^). 
Although many choices are possible, we take our multi-product formulas to be of the following form because they 
yield a large value of k while consisting of relatively few exponentials. 

Definition 1. Let k > and x ^ 1 integers, 7 a real number such that e'^^^'^^^ is an integer, and S^{X) a 
symmetric product formula approximation to U{X) obeying IIS'^(A) — C/(A)|| G 0{X^^^^). Then for any t M. we 
define the multi-product formula Mfc_^(<) as 

k+l 

Mu,^{t):=Y,CqS^{t/tqY^ (15) 



9=1 



wher 



^^-le^C^+i) ifq = k + l ^ ' 



r \ g2_e2^(fc+l) Ylj^q q2_j2 if Q < k 

{\Aj=i e27(fc+i)-j-2 ifq = k + l. 

We choose these values of £q because for sufficiently large 7 they guarantee that Ck+i is much larger in absolute 
value than all other coefhcients. This ensures a high success probability in Theorem 3 because k > jC'^+il/ |C'g| 
is large if |Cfc+i| exceeds the sum of all other |Cg|. 

The following lemma shows that M^^^ is a higher-order integrator than Sk- Quantitative error bounds are proven 
in the next section. 

Lemma 4. Let M^^^ be a multi-product formula constructed according to Definition 1. Then for A <C 1 we have 

|lMfe,^(A)-C/(A)|| eO(A^('=+x)+i). (^8) 

Proof. We follow the steps outlined in Chin's proof for the case where x = 1 [-2]. As shown in [19], a sufficient condition 
for a multi-product formula of the form Mk,^{X) = J2p CpS^iX/tp)^^ to satisfy ||t7(A) - Mk,^{X)\\ e 0(A2('=+x)+1) is 
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for C to satisfy the following matrix equation: 
/ 



1 



■2x 



-2x-2 



-2x 



*2 



/-2X 

«-2x-2 



„-2(fe+x-l) .-2(fc+x-l) „-2(fc+x-l) 
*-2 *3 



/-2x 

/-2x-2 



„-2(fc+x-l) , 







M 


C2 







C3 








\Ck+i J 



(19) 



This ensures that the sum of all Cq is 1 and the coefficients of all the error terms in the multi-product formula are 
zero up to 0(A^''''"*"'^-'~^). Denoting the matrix in the above equation by V, the vector C is then the first column of 
V~^. The matrix y is a generalized Vandermonde matrix, which can be explicitly inverted [2")]. The entries of C 
correspond to (13), which coincides with the values of Cq given in (17) for the values of £q in (16). Therefore the 
result of [ : ] shows that these values of Cq extrapolate an 0{X^^'^^) symmetric product formula into an 0(A^'^'''+'''+^) 
multi-product formula, as claimed. □ 

The following upper bound on the coefficients of the multi-product formula will be useful. 
Lemma 5. Ij e^'^^^'^^^ > 2k^, then for all 1 < q < k + 1, the coefficients Cq from Definition 1 satisfy 

\Cq\ < y2fc3/2e2Ml+logW/2-7) 

where 

A2 



V 



max 



Ae[o,i) (1 + A)i+^(1- A)i-^ 



0.3081. 



(20) 
(21) 



Proof. For any q < k + 1, Definition 1 gives 



Cq = 





_ e27(fe+i) 










q^ 


_ e27(fc+i) 








q' 



n 

je{l....,fc}\q 

n 



q2 _ p 



q 



{q + j){q-i) 



]e{i,...M\q 

2(k-i) 2q-q\ (-1)'--^ 
q2 _ e27(fe+l) + ky. {q - l)l{k - q)l 

(_l)fc-92q2fc+2 

" {k + q)l{k - qy.{q^ - e^-rik+i))- 
Using 6^'''^*^+^) > 2fc2 > 2q2, we have the bound 

4^2fe+2g-27(fc + l) 



We proceed by using the lower bound [20] 



{k + qy.{k-qy. 



i! > \/2^n"e-("+i/^3^ 



(22) 



(23) 



(24) 



We will use this bound differently to estimate \Cq\ for the cases where q < k and q = k. Using (24), we lower bound 
both factorial functions in (23) as follows: 



4„2fc+2p-27(fc+l)+2fe+2/13 
1(7 I < ^ 

27^^/fc(fc + g)fe+9(fc - q)k-q ' 



(25) 



Here we have used the fact that \/ (k + q){k — q) > \fk for g < fc — 1. Now introduce a parameter A such that q = fcA. 
We simplify (25) using this substitution and divide the numerator and denominator of (25) by k'^^ to find 



\C,\ < 



2^3/2g2fc(l-7)-27+2/13 



(1 + A)i+^(1-A)i-^ 



< 



2^3/2g2fc(l-7)-27+2/13 



(26) 
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We then simplify this expression and find that 



2^3/2g2fe(H-log(»,)/2-7) 



l^.l ^ ' (27) 

which for 7 > is bounded above by 

< A:3/2e2Ml+log(^)/2-7) < ^^3/2g2fc(l+log(^)/2-7)_ (28) 

Our bound for the case where g = fc is found using similar (but simpler) reasoning. We first substitute q = k 
into (23), use the lower bound in (24) to remove the factorial function, and simplify the result to find 

\Ck\ < 2fc3/2e2*^(l-l°s(2)-7)/^ < ^fc3/2g2fc(l-log(2)-7)^ (29) 

which is less than the value in (28) because -0.6931 « - log 2 < log(77)/2 sa -0.5886. Therefore (20) holds for all 
q < fc as required. □ 

The value of 7 used in Mk^^ '^^^ chosen to minimize the probability of a subtraction error. The following lemma 
relates the value of k to 7, allowing us to use Theorem 3 to find a value of 7 that ensures a sufficiently small probability 
of a subtraction error. 

Henceforth we assume that fc > 1. This is because k = corresponds to an ordinary product formula and the 
bounds that we prove for multi-product formulas can tightened by excluding this case, which is also already well 
analyzed [II, IS]. 

Lemma 6. Let M^^^ be a multi-product formula as in Definition 1 and let k be defined for M^^^ '^^ ™ Theorem 3. 
Then if2k^ < e^'^f'^+i) and k>0, we have 

K > 2-l/2g-2/c(l+log(,,)/2-7)-log(fc=/2) (3Q-) 

where rj is defined in (21). 

Proof. According to Theorem 3, we have k = S]+/I]_, where S+ is the sum of all Cq with positive coefficients and S_ 
is the absolute value of the corresponding negative sum. A lower bound on k is therefore found by dividing a lower 
bound on S_|_ by an upper bound on E_. 

Using the expression for Cq from Definition 1, we have 



g27(fe+l) ^ I 



Cfe+1 - Jl g27(fc+l) _ j2 ~ II X _ j2g-27(fc+l) ■ (^^) 



The denominators on the right hand side of (31) are positive under the assumption that 2k'^ < e'^'^^^^^\ which ensures 
that ;;2g-27(fc+i) ^ and simplifies the subsequent results of Corollary 7. We also have Ck+i > 1 because each 
denominator is less than 1. Since Ck+i > 0, we have > Ck+i > 1, and therefore 

«>^. (32) 

Next we provide an upper bound for S_. Since Ck+i > 0, we have 

fe 

S_<^|C,|. (33) 

q=l 

An upper bound for I]_ can then be obtained directly from upper bounds for max^^fc+i \Cq\. 
Using Lemma 5, we have 

k 

< ^ %/2 A;3/2e2Ml+log(r,)/2-7) < ^ ^5/2g2fc(l+log(^)/2-7) _ (34) 
9=1 

We substitute this inequality into (32) to obtain 

K > 2-l/2^-5/2g-2fc(l+log(.)/2-7) (35) 

as claimed. □ 
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In fact, the bound of Lemma 6 is nearly tight, in that k decays exponentially with fc if 7 < 1 + log{i])/2, as we 
discuss in more detail below. Thus the success probability of our algorithm decays exponentially if 7 is too small. 
Consequently, we will find that unlike the classical case, our quantum algorithm does not provide poly-logarithmic 
error scaling. 

The following corollary provides a sufficient value of 7 to ensure that the probability of our algorithm making a 
subtraction error is small. 

Corollary 7. Let Mk^x ^ multi-product formula as in Definition 1, let k be defined for M^^x '^■^ Theorem 3, and 
let S < 1. Furthermore, suppose k > and 



2fc \ S 



Then P_ <S and fc2e-27(fc+i) < 1/2. 



Proof. Without loss of generality, we can take k > 1 for our subtraction step because J^q Cq = \ and hence k, = 
E+/E_ > 1. This observation and the result of Theorem 3 imply that the probability of failing to perform the 
subtraction step in our implementation of Mfe.-^ satisfies 

P- < 4/k. (36) 

Eq. (36) and the bounds on n in Lemma 6 give P^ < S provided 

4y2e^''^^+'°*5(''^/^"''^+'°8i^'''^''^ < 6. (37) 
We obtain our sufficient value of 7 by solving (37) for 7, giving 

7>l + log(7y)/2 + ^log(^Miy (38) 

Eq. (38) also imphes that fc^e-^TC^+i) < 1/2. For < (5 < 1 and 7 saturating (38), it is easy to see that fc^g-^TC^+i) 
is a monotonically decreasing function of k, and therefore achieves its maximum value at fc = 1, the smallest possible 
value of k. We find that k^e'^'''-''+^^ < 1/2 at A: = 1 for (5 = 1, and therefore the condition fc2e-27('=+i) < 1/2 is 
automatically implied by our choice of 7 in (38). □ 

The value of 7 given by Corollary 7 is tight up to 0{k~^ logfc). This is illustrated in Figure 2, which shows k as a 
function of fc for 7 = 7c := 1 + log(?7)/2 and two slightly perturbed values of 7 centered around 7c. We see that small 
deviations away from 7 = 7c lead to either exponential growth of k or exponential convergence of k to 1. Thus our 
lower bound for 7 cannot be significantly improved. 



IV. ANALYSIS OF SIMULATION AND ERRORS IN MULTI-PRODUCT FORMULAS 



The results of the previous section show how k scales with the number of terms in the multi-product formula used 
in the simulation. We now expand on these results by bounding the approximation errors incurred by using the 
multi-product formula. We also present an error correction method to ensure that our implementation fails at most 
a constant fraction of the time. Our error bounds are established as follows. First, we estimate the error in multi- 
product formulas that utilize high-order Lie-Trotter-Suzuki formulas. Second, we discuss how these multi-product 
formulas are implemented. Finally, we estimate the inversion error for the resulting multi-product formulas and bound 
the average error resulting from a given step. 

Lemma 8. Let Mk,ki^) satisfy Definition 1 for evolution time A > 0, Zei |Cg| < 2 for all q ^ 1, . . . ,k + 1, and let 

\\UiX) - MkAm < {2m{b/it-^h\f^+\ (39) 

The proof of Lemma 8 requires upper bounds on the remainder terms of Taylor series expansions. Let R£(/) 
denote the remainder term of the Taylor series of a function / truncated at order I. The following lemma bounds the 
remainder term for an operator exponential. 
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FIG. 2: Scaling of k with k for three values of 7 centered around 7c 1 + log(T;)/2. The data show that k approaches 1 
polynomially quickly if 7 = 7^, whereas a slight increase in 7 causes k to grow exponentially and a slight decrease causes k to 
converge to 1 exponentially with k. 



Lemma 9. Let a, € K for j = 1, . . . , M and suppose \\Hj\\ < h. Then 



M 



Proof. Using the triangle inequality and sub-multiplicativity of the norm, we find 



M 



as claimed. 

Now we are ready to prove Lemma 8. 
Proof of Lemma 8. Lemma 4 implies that 



(Af / CO N 
n EH«,^,i)7p! 

(M / 00 \ 
X{[Y,{\a,\htY/p\ 

I I"" 

= Re exp ^ \aj\ht 



< 



E 

p=e+i 



i+1 



M 

■ exp I 



(40) 



(41) 



(42) 



□ 



||A4,fc(A)-C/(A)|| eO(A*'=+^), 



(43) 
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so the approximation error is entirely determined by the terms of order A^'^^^ in Mk.k and U . If we remove the terms 
in the Taylor series of Mk^k and U that cancel, the remainders, R4fc(Mfe^fe(A)) and R4fc(C/(A)), determine the error via 

||Mfc,fe(A) - U{\)\\ = ||R4fe(Mfc,fc(A)) - R4fe(C/(A))|| 

< ||R4fe(Mfc,fc(A))|| + ||R4fc(C/(A))||. (44) 



Lemma 9 implies that 

\\Ti.,ATT(\\\\\ < ' 

(4fc + 1)! 



R4.(C/(A))||< ^"'^^^''^'^-''^ 



'A 

<\-m{b/if-^h\] . (45) 



3 

4mfc(5/3)''- 



The second inequality in (45) follows from the assumption that h\ < j~]7[^^^k^: which implies that exp(m/iA)/(4fc + 
1)! < 23/4/5! < 1. 

The definition of Mk^^ implies 



fe+i 

||R4fc(A4,fc(A))|l < ^ \Cg\ \\lUk{Sk{X/eqY')\\ ■ (46) 

9=1 

Thus we upper bound ||R4/i;(S'/i;(A/p)P)||. This bound follows similar logic to the bound for U(X), but the calculation 
is slightly more complicated because 5*^ is the product of many exponentials. Specifically, 

Sk{X/p)= n e-'"^^'^'-'^/P, (47) 



=1 



where qk,i is the ratio between X/p and the duration of the £ exponential in Sk- Lemma 9 gives 

l|R4.(5.(A/p)^)|| < (2^5^- W./^A)4^+^ ^,^,.-. 

(4k + 1)! 

Using the upper bound g^/ < 2k/ 3'' from Appendix A of [2i)], we have 

IIR4.(^.(A/p)^)|| < (49) 



This bound can be simplified using (4fc + 1)! > ^4*^+^5! for fc > 1 (a consequence of (24)) and the hypothesis that 
^mk{5/3)''-^hX < log(2), giving 

4fe+l 



(50) 



\\Rik{Sk{X/pr)\\ < I Qm(5/3)'^-i/iA^ 
Using this result and the assumption that jC^j < 2 in (46) gives 

||R4fe(Mfe,,(A))|l < ^^^^ Qm(5/3)'=-i/iA) . (51) 



Combining (44), (45), and (51) gives 



||C/(A)-A4,fc(A)|| < ( 1+ ^ ' \ [-mlS>/if-^hX 



4fc+l 



< (2m(5/3)'=-i/iA)4'=+i, (52) 
proving the lemma. □ 
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A useful consequence of performing M^^k using a single subtraction step is that if a subtraction error occurs, the 
simulator performs the operation 



Ek{\): \'<P) ^ 



This error operation can be approximately corrected because, as Blanes et al. proved [ I'i, Theorem 1], 



(53) 



(54) 



Since the coefficients \Cq\ are all positive, Theorem 3 shows that the approximate correction operation Ek{—X) can be 
performed with success probability close to 1 provided A is small. The following lemma states that the error incurred 
by approximately correcting subtraction errors is at most equal to our upper bound for the approximation error for 
Mfc,fc(A). 



Lemma 10. Let -Efe(A) act as in (53), where Cq and £q are given in Definition 1. If 2mk{5/3f-'^hX < 1/2, then 

(a - EkhX)EkiX)) m < (2mA:(5/3)'=-i/iA)''+' . 



max 



(55) 



Proof By (53) and (54), 



max||(a-Sfc(-A)Sfc(A)) 



m£,xWRik+i{Eki-X)EkiX)\^))\\ 

IV') 



< 



R4fe+i (EpE, \Cp\\Cq\Sk{-x/prSk{x/qr 



(56) 



We then follow the same reasoning used in the proof of Lemma 9. By the triangle inequality, the norm of the remainder 
of a Taylor series is upper bounded by the sums of the norms of the individual terms in the remainder. This can be 
bounded by replacing the exponent of each exponential in Sk with its norm. We use similar reasoning to that used 
in (49) to find 



||R4fe+i {Ski-X/pfSkiX/qm < R4fe+i (^ei"'H5/')'-'hx^^ 
so the numerator of (56) satisfies 

(5]^|Cp||C,|5,.(-A/p)f5fc(A/g)M < ^5]|Cp||C,|R4fe+i 
\ p g / p q 



(57) 



|mfc(5/3)'="^/iA 



(|mfc(5/3)'=-i/iA)''+' et"'=(5/3)'=-^''^ 



(58) 



where C is a vector with entries Cp, so ||C||i — Ep IOjI- Our assumptions imply |mfc(5/3)'^ ^hX < 2/3 < log(2), so 

,2(|mfc(5/3)'=-i/iA)^^+' 



5]5]|Cp||C,|5fc(-A/p)''5fc(A/g)^ 



p q 



<\\crv 



(4fc + 2)! 



< ||C||22(f ™(5/3)'=-i/.A)^'=+^ 



(59) 



where the last inequality results from using Stirling's approximation as given in (24). 
The denominator of (56) can be lower bounded as follows: 



mm 

10) 



J2\Cp\Sk{x/epY 



mm 



^ |Cp|(e-»^* - (e-^* - Sk{X/£,Y'')M 



> \\C\\i l-max|le-^^*-5fe(A/^p)^-| 



-iHt 



SkiX)\\). 



(60) 
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Since 2mk{b/?,f-^h\ < 1/2 < 3/(472), Theorem 3 of [20] implies that Ije"'^* - S'fc(A)|| < 2{2mk{5/3)''-^hX)'^''+^ . 
Using 2mfc(5/3)''-i/iA < 1/2 and fc > 1 we have that 



-iHt 



e 



Sk{X)\\ < I, (61) 



implying 



mm 



>l\\Ch. (62) 



Combining this with our upper bound on the numerator gives 



< (2m(5/3)'=-i/iA)4^+2 (63) 

as claimed. □ 

We simulate U{t) using r iterations of Mk^k{t/r) for some sufficiently large r. Our next step is to combine Lemma 8 
and Lemma 10 to find upper bounds on r such that U{t) is approximated to within some fixed error. We take S = 1/2, 
i.e., we accept a maximum failure probability of 1/2 for each multi-product formula. We then sum the cumulative 
errors and use the Chernoff bound to show that, with high probability, the simulation error is at most e. These results 
are summarized in the following lemma. 

Lemma 11. Let Mj^ be a multi-product formula given by Definition 1 with \Cq\ < 2 for all q < k + 1. Let 7 be 
chosen as in Corollary 7 with (5 = 1/2 and let the integer r satisfy 

(2^(5/3)^2^^ 

(e/5)i/4fe ^"^^ 

for e < mhtk^'^^ . Then a quantum computer can approximately implement U{t) as Mk,k{t/rY with error at most e 
and with probability at least 1 — e~^^^^ , assuming that no addition errors occur during the simulation, while utilizing 
no more than 5r subtraction attempts and approximate inversions. 

Proof. First we bound the probability of successfully performing the subtraction steps given a fixed maximum number 
of attempts. We simplify our analysis by assuming that the simulation uses exactly 3r subtractions, corresponding 
to the worst-case scenario in which 2r inversions are used. We want to find the probability that a randomly chosen 
sequence of subtractions contains at least r successes, correctly implementing the multi-product formula. The proba- 
bility that a sequence is unsuccessful is exponentially small in r because for i5 = 1/2, the mean number of failures is 
IJ, = 3r/2, which is substantially smaller than our tolerance of 2r failures. By the Chernoff bound, the probability of 
having more than 2r failures satisfies 

Pt{X > 2r) < e-^((i+")i°g(i+")^") < e""^/", (65) 

where 1 + a = 2r / fi — 4/3. 

If we attempt the subtraction steps in our protocol 3r times and fail 2r times, then 5r subtractions and approximate 
inversions must be performed, because every failure requires an approximate inversion. We bound the resulting 
error using Lemma 8, Lemma 10, and the subadditivity of errors. These lemmas apply because the requirement 
^mk{5/3Y~^ ht /r < log(2) is implied by our choice of r and the assumption e < mhtk~^''. By this argument, the 
simulation error satisfies 

\\Mk.kit/rY - Uit)\\ < 5r{2m{5/3Y-^ht/rY''+^ (66) 

with probability at least 1 — e^'"/^'^, where Mk^kit/r) denotes the operation performed by our non-deterministic 
algorithm for the multi-product formula Mk^k- The assumption e < mhtk~'^^ and the value of r from (64) imply that 
\\Mk,k{t/ry - U{t)\\ < e as required. □ 

Lemma 11 assumes an error tolerance of at most mhtk~^^ , which may appear to be very small. However, our 
ultimate simulation scheme has k G 0{^y\og{mht/e), so in fact the error tolerance is modest. 

Now we are ready to prove a key lemma that provides bounds on the number of exponentials used by the simulation 
in the realistic scenario where both addition and subtraction errors may occur. Our main result. Theorem 1, follows 
as a simple consequence. 
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Lemma 12. Let M^^k be a multi-product formula for H = Yl^=i^j '^^ ^'^ Definition 1, with k > 1. Let M^^k be 
the implementation of Mk^k described above. Let e be a desired error tolerance and let (3 he a desired upper bound 
on the failure probability of the algorithm. Then there is a simulation of U{t) = e^'^* that has error at most e with 
probability at least 1 — /3 using 

TVexp < 1000TO5'=-iA:^/4e(i+'°s('?)/2)V (67) 

exponentials of the form e~*^^*, where 

l-l^l logrexp([l + log(r;)/2 + ^ \og{2{2kfl^)]k)-\ , 
2. e^Ta\n{l,e,l3,mhtk~'^^), 



3. r ■ 



f^"'^%v"i""^^ l3log(2/^)}^ 



Proof. We approximate U{t) as a product of r multi-product formulas, each implemented using the operation Mk^k- 
We can suppose that the sequence of r multi-product formulas is implemented using at most 5r subtraction and 
inversion steps according to Lemma 11. As fc unitary operations are combined in each step, we must implement a 
total of brk unitary operations. Definition 1 implies that each of these unitaries is composed of at most e^^^'^^'^ Suzuki 
integrators Uq. Each Uq is a product of 2r7i5'^~^ exponentials of elements from {Hj}. Thus, if the algorithm succeeds 
after performing this maximum number of subtractions and inversions, we have 

A^cxp < 10m5'="ifce'^('=+iV 

< 10m5'=-lfc(2"/4e2fc5/4g(i+log(^)/2)fc ^ 

< 1000m5'=-ifc9/4e(i+'°s('))/2)fe^ (68) 

where we have used 7 < 2 and e'^'' < 2(2'^/'*fc^/'*e'=(i+^°s('')/2)). Equation (67) then follows by substituting the value of 
r assumed by the lemma, which guarantees that the simulation error is less than e when the simulation is successful 
because it exceeds the value of r from Lemma 11. We choose r to be larger because it simplifies our results and 
guarantees that the probability of an addition error is at most e/2. 

Lemma 11 requires Cq < 2, which we have not explicitly assumed. Substituting the above value of 7 into the upper 
bound for E_ in (34) shows that X]g=i — 1' l^ql — 1 &\\ q = 1, . . . , k. Our multi-product formula satisfies 
J2„ Cq — I, so our choice of 7 ensures Ck+i < 2. Thus Cq < 2 for all q. 

Lemma 11 implies that the probability of the simulation failing due to too many subtraction errors is at most 
g-r/13 However, it does not address the possibility of the algorithm failing due to addition errors. There are at most 
5r addition steps, so by Theorem 3 and the union bound, the probability of an addition error is at most 

5rP+ < ——. (69) 



By the definition of A in Theorem 3, 



A = 



max 
q^q 



.xR2fc {Skit/kqrp - Sk{tlkq,rf^') 



(70) 



Using Lemma 9 (similarly as in (49)), the triangle inequality, and |m/i:(5/3)^ ^ht/r < log(2), we have 

(2fc+l)!" 



(|mA:(5/3)^-i/it/r)2fe+i ^ ^ 

A < 4^^ —-^ — . (71) 



Substituting the assumed value of r into (71) gives 



20fc2fe+i / i 
A < T^TT^^TTT^ .^../o^.-1.. ■ (72) 



(2/c + l)!32fe+i V4to(5/3)'=-iM 
Substituting this bound into (69) and using (24), we find that the total probability of an addition error is at most 



((2fc-hl)!)234'=+2 \Amk{b/?,Y-^htJ " 67r(6/e)4fc+2e-2/i3 \^4mA;(5/3)'=-i/ii 
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Since e < mhtk , this implies 



20? 



67r(6/e)4fc+2e-2/i3 \^4/c(5/3)'^-i 



l/4fe 



e e 
130 2' 



(74) 



Here the second inequaUty follows from the first by substituting k = 1, since the middle expression is a monotonically 
decreasing function of k. 

The total probability of success Pg satisfies 



> 1 - - - 
2 

Using r > 13 log(2//3) and e < /?, we find > 1 — /3 as claimed. 



(75) 

□ 



As in [11], there is a tradeoff between the exponential improvement in the accuracy of the formula and the expo- 
nential growth of Mfe fe with k. To see this, note that apart from terms that are bounded above by a constant function 
of fc, -/Vcxp is the product of two terms: 



A^oxp € O 



G O 



^2^9/4g(l+log(r,)/2+log(25/3))/c 



[mht/e 



~il/4fe 



m'^k^^^e 



4„2.54fc 





mht 


l/4fe\ 




min(e, /3) 





(76) 
(77) 



Here we have not included the 0(log(l//3)) term from r because log(l//3) e 0(l/^)i/''^ For comparison, the results 
of [11] and [20] have the following complexities: 

iVe.pGO([m2e3-22'=][m/iVe] (78) 

iVexp e O ([m^/fce^-^fc] [m/ii/e]^/''=) , (79) 

respectively. The tradeoff between accuracy and complexity as a function of k is more favorable in our setting than 
in either of these approaches. Finally, we give a detailed analysis of the tradeoff. 

Proof of Theorem 1. Neglecting polynomially large contributions in (76), we see that the dominant part of (76) is 

g(l+log(i,)/2+log(25/3))fc+4^ log(Tn/it/e)_ ^gQ-j 

It is natural to choose k to minimize (80). The minimum is achieved by taking k — k^pt, where 

0.3142 v/log(m/ii/?). (81) 



^opt — 



1 / \og{mht/i) 
2V l + log(??)/2 + log(25/3) 



Using this fc, we find that 



iVexp e O {k%tm^hte'-^V^°s(mht/e)^ 



(82) 



We have e = min(l, e, /3, m/iifc ^'^), so e depends implicitly on fc. However, we now show that this term can be 
neglected in the limit of large mht/e. In this limit, we have 



fcopt ^ ^ ( \og{mht/e) V I ^ ^^^^^^^^ 



mht 



mht 



(83) 



The above follows since lima;^oo log(a;)'^V^^(^/a; — Q for any c > 0. Correspondingly, mhtk^p^""^ € ^{^)- We therefore 
conclude that this term can be neglected asymptotically and that e S ri(e), so we can replace e with e asymptotically. 
The result then follows by substituting fcopt into (82) and dropping all poly-logarithmic factors. □ 

The parameters e and /3 can be decreased to improve the simulation fidelity and success probability, respectively. 
The cost of such an improvement is relatively low, although it is not poly- logarithmic in 1/e and 1//3. If the initial 
state of a simulation can be cheaply prepared and the result of the computation can be easily checked, it may be 
preferable to use large values of e and (3 and repeat the simulation an appropriate number of times. Then the Chernoff 
bound implies that a logarithmic number of iterations is sufficient to achieve success with high probability. However, 
this may not be possible if the simulation is used as a subroutine in a larger algorithm (e.g., as in [(>]). 
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V. CONCLUSIONS 



We have presented a new approach to quantum simulation that implements Hamiltonian dynamics using linear 
combinations, rather than products, of unitary operators. The resulting simulation gives better scaling with the 
simulation error e than any previously known algorithm and scales more favorably with all parameters than simulation 
methods based on product formulas. Aside from the quantitative improvement to simulation accuracy, this work 
provides a new way to address the errors that occur in quantum algorithms. Specifically, our work shows that 
approximation errors can be reduced by coherently averaging the results of different approximations. It is common 
to perform such averages in classical numerical analysis, and we hope that the techniques presented here may have 
applications beyond quantum simulation. 

It remains an open problem to further improve the performance of quantum simulation as a function of the error 
tolerance e. In particular, we would like to determine whether there is a simulation of n-qubit Hamiltonians with 
complexity poly(n, log ^). Classical simulation algorithms based on multi-product formulas achieve scaling polynomial 
in log^, but they are necessarily inefficient as a function of the system size. Our algorithms fail to provide such 
favorable scaling in ^ due to the sign problem discussed in Section III. A possible resolution to this problem could 
be attained by finding multi-product formulas with only positive coefficients and backward timesteps. Such formulas 
are not forbidden by Sheng's Theorem [--j] since that result only applies to multi-product formulas that are restricted 
to use forward timesteps [18]. New approximation-building methods might use backward timesteps to give multi- 
product formulas that are easier to implement with our techniques. Conversely, a proof that Hamiltonian simulation 
with poly (n, log-) elementary operations is impossible would also show that Lie-Trotter~Suzuki and multi-product 
formulas with positive coefficients cannot be constructed with polynomially many exponentials, answering an open 
question in numerical analysis. 

Finally, we have focused on the case of time-independent Hamiltonian evolution. One can consider multi-product 
formulas that are adapted to handle time-dependent evolution, as discussed in [_ , ]. It is nontrivial to use such formulas 
to generalize our results to the time-dependent case because of difficulties that arise when the Hamiltonian is not a 
sufficiently smooth function of time. Further investigation of these issues could lead to developments in numerical 
analysis as well as quantum computing. 
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Appendix A: Optimality of the linear combination procedure 



The goal of this appendix is to show that no protocol for implementing linear combinations of unitary operations 
in a large family of such protocols can have failure probability less than 

(Al) 



where k is defined in Theorem 3. Specifically, we consider protocols of the form shown in Figure 3. Our result shows 
that the protocol of Theorem 3 is optimal among such all protocols in the limit of small A (i.e., when the unitary 
operations being combined are all similar). 

Theorem 13. Any protocol for implementing V — X]g=o ^q^q using a circuit of the form of Figure 3 must fail with 
probability at least 4k/{k + 1)^. 



Proof. For convenience, we take k to be an integer power of 2. We can generalize the subsequent analysis to address 
the case where fc -I- 1 is not a power of 2 by replacing k + 1 hy k' + I = 2r'°S2('^+i)l and taking Cp = for p > k. 
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|0) 

|o) H 
|o) 



A 



B 



FIG. 3: A general circuit for implementing a linear combination of fc + 1 unitary operators using unitary operations A and B. 
We assume for simplicity that + 1 is an integer power of 2. This circuit corresponds to preparing the ancilla states in an 
arbitrary state (specified by A) and measuring them in an arbitrary basis (specified by B). 



Observe that the circuit in Figure 3 acts as follows: 



m— 

^ y^^Ajn,Q\m)Um\lp) 

m 

I-*- y^^Bn.mAm,o\n)Um\lp). 



(A2) 



Furthermore, we can modify B to include a permutation such that desired transformation occurs when the first 
register is measured to be zero. (Orthogonality prevents us from having more than one successful outcome, as we 
show below.) The implementation is successful if there exists a constant K > Q such that 



KC„ 



(A3) 



for all m. 

Our goal is to maximize the success probability, which is equivalent to maximizing K over all choices of the unitary 
operations A and B satisfying (A3). We can drop the implicit normalization of the matrix elements of B and A by 
defining coefhcients &o,m and a^.o such that 



0,m 



•— Am,0- 



Using these variables, we have 

\KC^\ = 



I^O.mOm.ol 



(EJa.,oP)(EJV 



< 



\bo,mO"m,o\ 



\C„ 



(A4) 
(A5) 

(A6) 



where the bound follows from the Cauchy-Schwarz inequality. This bound is tight because it can be saturated by 
taking am,o — ^o,m = VCm- The probability of successfully implementing the multi-product formula is 



l-P^ = \\Y,KC,U,\^) 



< 



' E,-fi 
E.IQI 



where and I]_ are defined in Lemma 6. We have k := so 

2 



1 - P_ < 



K-1 
K+1 



(A7) 



(A8) 
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which imphes that the failure probabiHty satisfies P_ > 4k/(k + 1)^ as claimed. 

It remains to see why it suffices to consider a single successful measurement outcome. In principle, we could imagine 
that many different measurement outcomes lead to a successful implementation of V. Assume that there exists a 
measurement outcome v ^ 0'°S2 ^ such that the protocol also gives the same multi-product formula on outcome v. 
If both outcomes are successful then, up to a constant multiplicative factor, the coefficients of each Uq must be the 
same. This occurs if there exists a constant F 7^ such that for all q, 

Bo.qA-qfi — TBy^qAq^o, (A9) 

i.e., if Bo q = TBy q (note that Aq^ must be nonzero provided Cq ^ 0). This is impossible because B is unitary and 
hence its columns are orthonormal. Consequently, we cannot obtain the same multi-product formula from different 
measurement outcomes. □ 

The above proof implicitly specifies an optimal protocol for implementing linear combinations of unitaries. However, 
we do not use this protocol in Theorem 3 because it is difficult to perform a correction if the implementation fails. 
When the method of Lemma 2 fails to implement a difference of nearby unitaries, the desired correction operation is 
a sum of nearby unitaries, so it can be implemented nearly detcrministically. The correction operation may not have 
such a form if we use the protocol implicit in the above proof. 

One simple generalization of the form of the protocol shown in Figm-e 3 is to enlarge the ancilla register and allow 
each unitary in the linear combination to be performed conditioned on a higher-dimensional subspace of the ancilla 
states. It can be shown that a certain class of protocols of this form also do not improve the success probability. 
Whether protocols of another form could achieve a higher probability of success remains an open question. 
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